
double dist=0.0;
double ibin=0.0;
int upb=0;
int lowb=0;
int iphy=0;

for (int i = 0; i < *nfocal; i++){
    for (int j= 0; j < *ntotal; j++){
    	if( (fsp[i]==sp[j]) ||
    	    (fx[i]>(x[j]+*rmax)) || (fx[i]<(x[j]-*rmax)) ||
    	    (fy[i]>(y[j]+*rmax)) || (fy[i]<(y[j]-*rmax)) ||
    	    (x[j]==fx[i] && y[j]==fy[i] )){
    	    continue;
    	}
    	dist=pow(pow(x[j]-fx[i],2)+pow(y[j]-fy[i],2),0.5);
    	if(dist > *rmax ){
    	    continue;
    	}
    	ibin=dist / *step;
    	upb=ceil(ibin)-1;
    	lowb=floor(ibin);
    	if(lowb!=0)
    	    lowb--;
    	 bincout[upb]++;
    	 bincout[lowb]++;

        iphy = sp[j]*(*nsp)+fsp[i];
    	binsum[upb]+=phyd[iphy];
    	binsum[lowb]+=phyd[iphy];
    }
}

